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Abstract 

The  behavior  of  the  liquid-solid  interface  during 
resolidification  of  a  condensed  coppei — titanium  crystal  is 
studied  by  atomistic  simulations.  A  system  of  1080  atoms 
(540  of  each  type)  was  simulated.  The  interface  was  created 
using  monte  carlo  techniques  and  resolidification  was 
performed  by  molecular  dynamics.  Resolidification  was 
attempted  at  temperatures  of  1200  K,  1100  K,  1000  K,  900  K, 
818  K  and  773  K.  The  liquid-solid  interface  was  found  to 
move  only  at  773  K,  resulting  in  a  chemically  disordered 
crystal.  The  interface  velocity  observed  in  the  simulation 
is  6.6  m/s  which  agrees  well  with  the  theoretical  value  of 
5.5  m/s.  The  low  interface  velocity  indicates  the 
recrystallization  is  governed  by  diffusion-limited  growth. 


vi 


Atomistic  simulation  is  one  way  to  predict  the  results 


of  various  resolidification  processes.  In  this  study,  a 
computer  simulation  is  used  to  investigate  the  behavior  of  a 
liquid-solid  interface  during  rapid  resolidification  of  the 
i ntermeta 1 1 i c  compound  CuT i .  An  interface  is  created  in  a 
simulated  CuT i  crystal  and  then  the  system  is  held  at 
specific  temperatures  for  several  thousand  timesteps  of 
0.002  picoseconds  each.  The  results  from  the  computer 
simulations  are  examined  to  see  if  the  system  has 
resolidified  and  the  interface  moved,  and  if  possible, 
calculate  the  interface  velocity.  The  crystal  may  grow  back 
chemically  disordered,  meaning  that  atoms  occupy  lattice 
sites  but  are  not  the  same  type  which  originally  inhabited 
the  location.  The  system  could  grow  back  amorphous  or  even 
undergo  homogeneous  freezing,  forming  a  solid  without  the 
interface  moving. 


Although  this  appears  to  be  a  straightforward  problem, 
all  simulations  and  experiments  performed  on  liquid-solid 
metal  interfaces  to  date  deal  with  single  metal  systems 
(Landman,  et  al,  1988a,  1988b).  The  simulation  of  these 
interfaces  in  alloys  is  a  new  region  of  study,  with  no 
previous  work  available  for  background.  Therefore  the 
intent  of  this  paper  is  not  to  arrive  at  a  specific 
conclusion  about  the  interface  velocity  (such  as  determining 
the  dependence  of  the  velocity  on  temperature)  nor  to 
develop  new  methods  of  analysis  or  simulation,  but  rather  to 
perform  initial  simulations  of  alloys  and  create  a  basis  for 
future  simulations. 

PreYio.y.s 

Several  atomistic  simulation  studies  on  the 
recrysta 1 1 i zat ion  of  metals  have  been  accomplished,  either 
using  simulations  only  (Landman  et  al.,  1988),  or  in 
coniunction  with  experimental  studies  (MacDonald  et  al., 
1985).  All  the  work  so  far  has  been  concentrated  on  single 
metal  crystals  but  the  results  from  those  studies  are  used 
as  a  yardstick  or  the  results  of  this  study. 

MacDonald,  et  al . ,  studied  the  crystallization  of 
silicon,  copper  and  gold  using  both  experimental  and 
computational  methods.  They  measured  interface  velocities 
of  up  to  100  m/s,  which  supports,  the  hypothesis  pure  metal 
recrystallization  is  governed  by  col  1 ision- 1 imited  growth. 
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Simply  put,  a  crystal  in  the  co 1 1 i s i on- 1 imi ted  growth  regime 
grows  quickly  since  any  atom  moving  into  a  lattice  site  can 
fill  the  site  with  no  adverse  effects  on  the  system  energy. 
The  maximum  interface  velocity  for  co 1 1 i s i on- 1 imi ted  growth 
is  the  speed  of  sound  in  the  material 

The  other  growth  regime  described  by  MacDonald  is 
d i f f us i on- 1 i m i ted  growth,  which  requires  a  more  specific 
r ear r angement  of  the  atoms  than  co 1 1 i s i on- 1 imi ted  growth. 

D i f f us i on- 1 im i ted  growth  is  prevalent  in  concentrated 
metallic  alloys  such  as  CuT i .  As  the  system  is  supercooled 
the  atoms  will  diffuse  into  lattice  sites,  thus  restricting 
r ecr y sta 1 1 i zat i on  velocity  by  the  diffusion  coefficient,  D- 

vD~\/\2/ Dl- Dl/\  (1) 

A.  is  the  lattice  constant  (-  one  interatomic  distance).  The 
diffusion  coefficient  for  liquids  is 


Dl 


MSD 

2nbt 


(2) 


MSD  is  the  mean  square  displacement  of  the  atoms  (see 
Chapter  l l l  for  an  explanation  of  MSD),  n  is  the  number  of 
degrees  of  freedom  of  the  system  and  At  is  the  t’me  interval 
over  which  the  MSD  is  calculated.  For  three  dimensional 
systems,  n  =  3. 

The  velocities  will  initially  increase  with  decreasing 
temperature  because  at  higher  temperatures  the  atoms  have 
enough  energy  to  break  out  of  lattice  sites.  At  some 
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temperature,  the  velocity  will  begin  to  decrease  with 
decreasing  temperatures  because  the  atoms  are  moving  slower 
and  take  longer  to  diffuse  into  the  lattice  sites.  Since 

o 

D  L~ \0~scm2/s  and  X-2  A  for  metals,  the  growth  velocities  in 
the  d i f f us i on- 1 imi ted  regime  are  on  the  order  of  5  m/s 
(MacDonald  et  al . ,  1985). 

Richardson  and  Clancy  performed  a  simulation  of 
resolidification  of  gold  and  copper  (Richardson  and  Clancy, 
1990).  They  simulated  laser  processing  of  metals  by  setting 
up  a  crysta 1 -vapor  interface  and  introducing  an  energy 
"pulse"  onto  the  crystal  surface  which  is  at  room 
temperature.  The  fluence  and  duration  of  the  pulse  (5 
mJ/cnrr 2  and  20  ps)  were  based  on  the  laser  pulses  used  in 
the  experiments  by  MacDonald,  et  al.  (see  previous  chapter). 
After  the  "pulse"  is  switched  off,  Richardson  and  Clancy 
calculated  the  interface  position  and  temperatures  as 
functions  of  time. 

Richardson  and  Clancy  calculated  interface  velocities 
at  three  temperatures :  0.22  Tm  (melting  temper atur e ) ,  0.53 
Tm  and  0.92  Tm .  Their  results  for  copper  yielded  interface 
velocities  of  82,  57  and  9  m/s  respectively.  They  also 
observed  some  disorder  in  the  recrystal  1 i zed  region  at  the 
highest  temperature  and  contraction  between  the  outermost 
crystalline  layers  after  resolidification.  The  results  they 
obtained  for  gold  did  not  agree  well  with  experimental 
values  obtained  by  MacDonald  et  al,  finding  the 
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resolidification  velocity  was  -  40  m/s.  This  discrepancy 
was  probably  due  to  deficiencies  in  their  potential  model 
(Richardson  and  Clancy,  1985:  335). 

Landman  et  al .  ran  two  simulations  on  silicon  to 
determine  the  effect  of  growth  direction  on  the  velocity. 

The  first  (Landman  et  al . ,  1988a)  simulated  crystal  growth 
on  a  (100)  interface.  They  found  that  crystallization  did 
not  occur  until  the  system  was  undercooled  by  150  K  and  then 
the  crystal  growth  proceeded  at  an  average  velocity  of  -  18 
m/s.  The  second  simulation  was  performed  on  a  (111) 
interface  (Landman  et  al . ,  1988b).  Again,  150  K 
undercooling  was  required  before  recrysta 1 1 i zat i on  was 
observed.  They  employed  two  different  cooling  rates  to 
their  system.  The  "fast"  rate  simulations  gave  a  growth 
velocity  of  -  14  m/s  and  the  "slow"  simulations  resulted  in 
a  growth  velocity  of  9  m/s.  The  flat  solid-liquid  surface 
was  maintained  during  the  growth.  The  crystallization  was 
discontinuous  on  the  (111)  surface;  the  crystal  would 
experience  periods  of  fast  growth  disrupted  by  intervals  in 
which  the  interface  would  not  move.  In  contrast,  the  (100) 
surface  grew  uninterrupted. 

LCDR  Charles  L.  Wood,  USN,  performed  atomistic 
simulation  of  the  melting  point  of  CuT i  as  part  of  a 
research  apprenticeship  at  the  Air  Force  Institute  of 
Technology.  He  used  the  same  computer  code  and  interatomic 
potentials  as  in  this  study.  Wood  found  the  simulation 
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yielded  a  melting  point  for  CuTi  of  1287  K  which  compared 
well  with  tne  experimental  value  of  1257  K.  Although  he  did 
not  study  interfaces,  his  results  indicate  the  potential 
functions  used  in  this  simulation  are  accurate. 

The  rationale  for  studying  interfaces  in  metals  may  not 
be  readily  apparent.  In  laser  processing  of  materials,  the 
laser  imparts  short  (pico-  or  nanosecond)  pulses  of  energy 
into  the  metal,  creating  large  temperature  gradients.  The 
surface  of  the  metal  melts  and  the  liquid-solid  interface 
moves  deeper  into  the  metal.  When  the  pulse  ends  the  heat 
flows  into  the  solid  substrate  causing  rapid  cooling  of  the 
melt.  The  melted  region  becomes  undercooled  below  the 
equilibrium  freezing  temperature,  causing  the  interface  to 
reverse  and  progress  toward  the  surface  of  the  metal .  The 
kinetics  of  the  heat  flow  determine  the  structure  of  the 
resulting  solid,  either  crystal  or  amorphous  (Chokappa  et 
al.,  1989).  Simulations  like  this  one  will  become  more 
important  as  new  material  processing  methods  are  developed. 

Out  1 ine 

Chapter  I  I  gives  a  background  of  atomistic  simulation. 
The  information  is  sufficient  to  give  a  basic  understanding 
of  the  mechanisms  of  atomistic  simulation.  The  specific 
computer  code  used  in  this  study  is  also  explained,  and  a 
brief  description  of  the  crystal  simulated  is  given. 

Chapter  III  describes  in  some  detail  the  procedure  used  to 
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simulate  resolidification  of  the  CuTi .  Chapter  IV  analyzes 
the  results  of  the  simulation  at  773  K  which  yielded  a 
recrystallizes  system.  Chapter  V  summarizes  the  work  and 
states  the  conclusions  drawn  from  the  simulation. 
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Method 


Atomistic  Simulation  (Sabochick,  undated) 

Atomistic  simulation  is  a  method  of  performing  computer 
experiments  on  materials  based  on  the  premise  that  if 
interactions  between  individual  atoms  are  known,  the  motions 
of  the  atoms  can  be  calculated,  thereby  simulating  small 
systems.  Atomistic  simulations  fall  into  three  categories: 
molecular  statics  (MS),  monte  carlo  (MC)  and  molecular 
dynamics  (MD).  Molecular  statics  is  used  to  find  stable 
energy  configurations  of  atoms  in  crystals,  and  is  not  very 
useful  in  non-equilibrium  studies  such  as  this.  Molecular 
dynamics  simulates  the  behavior  of  the  atoms  by  solving 
Newton's  equations  of  motion  for  the  system,  giving  the 
position  and  velocity  of  each  atom  in  the  system. 

For  atomistic  simulation  to  work  properly,  the 
interatomic  potential,  which  determines  how  the  atoms 
interact  with  one  another,  must  be  known.  The  interatomic 
potential  is  used  to  calculate  the  total  energy  of  the 
system.  Oh  and  Johnson  developed  an  approach  to  developing 
interatomic  potentials  based  on  the  embedded-atom  method 
( EAM )  which  is  used  in  this  study  (Oh  and  Johnson,  1988). 

In  EAM,  the  potential  energy  for  a  system  is  given  by 

UM-  I  (3) 

t  '  /<( 
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U tet  is  the  potential  energy,  r,  is  the  position  of  atom  i  , 


Ft Cp<)  is  the  energy  needed  to  imbed  an  atom  i  into  the 
electron  density  p,  and  4 >(/(r(y)  is  the  repulsive  potential 
between  atom  i  and  atom  j  separated  by  a  distance  rc/.  The 
subscripts  k  and  m  denote  the  species  of  atoms  i  and  j 
respectively.  The  electron  density,  p£  is  found  by 

Pco-X/^)  (4) 

!+<■ 

/,(r(/)  is  the  contribution  by  atom  j  to  the  electron  density 

to  the  electron  density  at  atom  i  at  a  distance,  rt)  (Oh  and 
Johnson,  1988:471-472).  For  the  potentials  used  in  the 
present  work,  the  chosen  form  of  the  repulsive  potential  is 

4)*,n(r </)"4>e  fcm(r t//r e>(tin)"Vfcm  (5) 

and  the  electron  density,  f,  is 

/m(r<y)-/..mexp(-pmrty/r.>mfn)  (6) 

where  and  v  are  constants  and  r,i(nm,  r.,*m  are  scaling  factors 
to  make  the  exponent  dimensionless  (Sabochick  and  Lam, 
1990:7-8) . 

The  force  on  the  atom  is  calculated  as  the  negative 
value  of  the  partial  derivative  of  the  total  potential 
energy  on  the  atom 


-  d£/(r() 


(7) 
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To  avoid  wasting  computer  time  by  calculating  the  forces 
between  atoms  which  do  not  interact,  atomistic  simulations 
set  up  a  'neighbor  list’  for  each  atom.  The  neighbor  list 
specifies  which  atoms  will  interact  with  each  other.  The 
list  is  periodically  updated  during  a  run  to  account  for 
atoms  moving.  Another  computer  time  saving  method  is  to  set 
the  force  on  atom  j  from  atom  i  equal  to  the  negative  force 
on  atom  i  from  atom  j.  This  avoids  calculating  the  same 
force  twice.  Once  the  force  between  the  atoms  are  known, 
the  equations  of  motion  are  solved.  The  equations  of  motion 
cannot  be  solved  analytically  for  realistic  interatomic 
potentials,  and  thus  are  solved  numerically. 

Because  the  forces  on  an  atom  depend  on  the  number  of 
neighbors  it  has,  atoms  on  a  vacuum  only  see  half  the  forces 
that  an  atom  on  the  interior  would  experience.  Since  it 
would  use  an  inordinate  amount  of  computer  time  to  run  a 
system  of  a  million  atoms,  boundary  effects  are  minimized  by 
the  use  of  periodic  boundaries,  in  which  the  system  is 
surrounded  by  copies  of  itself.  Thus,  atoms  on  the  edge 
’see’  an  atom  across  a  boundary  and  react  accordingly.  If 
the  atom  crosses  a  boundary  during  the  simulation,  its 
’image’  will  enter  the  simulation  volume  across  the  boundary 
opposite  the  one  crossed  by  the  initial  atom. 

Five  state  variables  are  present  in  each  molecular 
dynamics  simulation:  volume,  pressure,  number  of  particles, 
temperature  and  total  energy.  For  any  MD  simulation,  three 
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of  these  variables  are  held  constant.  In  this  study,  the 
volume,  number  of  atoms  and  system  temperature  are  fixed, 
and  the  other  two  properties  are  calculated.  The  potential 
energy  is  calculated  using  equation  (1).  The  kinetic  energy 
(K.E.)  and  temperature  are  related  by 

K.E.- |  k*N*T  (8) 

where  k  is  Boltzmann's  constant,  N  is  the  number  of  atoms  in 
the  system  and  T  is  the  absolute  temperature.  The 
temperature  is  held  constant  using  the  extended  system 
method  (Nose,  1984). 

After  the  forces  and  temperature  are  determined,  the 
pressure  of  the  system  is  calculated 

NkT  1  ydu(  r„l 

V  31/ ,4}  dr  i, 

in  this  equation,  V  is  the  volume  of  the  system,  the 
derivative  of  the  potential  energy  is  the  magnitude  of  the 
force  of  atom  j  on  atom  i  and  rt/  is  the  distance  between 
atoms  i  and  j.  If  the  pressure  is  held  constant  rather  than 
the  volume,  the  system  will  expand  or  contract  with  an 
increase  or  decrease,  respectively,  in  pressure. 

Molecular  statics  calculates  the  forces  between  atoms 
in  the  same  manner  as  MD,  but  instead  of  solving  the 
equations  of  motion,  MS  uses  energy  minimization  routines. 

At  the  end  of  a  molecular  statics  run  the  atomic  positions 
should  be  such  that  the  potential  energy  of  the  system  is 


minimized.  Simply  put,  as  an  atom  passes  through  its 
position  of  minimum  potential,  its  acceleration  and  velocity 
of  opposite  sign.  When  the  product  of  the  velocity  and 
acceleration  becomes  negative,  the  velocity  is  set  to  zero. 
The  atom  begins  to  regress  along  its  path  until  the  product 
of  the  acceleration  and  velocity  becomes  negative  again,  and 
so  on.  Since  discrete  time  steps  are  used,  the  chance  of 
the  atom  coming  to  rest  precisely  at  its  lowest  potential 
configuration  is  low.  Hence,  a  molecular  statics  run  is 
stopped  when  the  maximum  force  on  any  atom  is  below  a  given 
to  1 erance . 

Monte  carlo  methods  are  usually  used  for  a  system  of 
atoms  in  a  fixed  volume  at  a  constant  temperature.  The 
Metropolis  monte  carlo  technique  goes  as  follows.  The  total 
energy  of  the  system  is  calculated  using  the  interatomic 
potentials.  An  atom  is  chosen  at  random  and  moved  to  a  new 
position.  If  the  energy  of  the  system  decreases  the  move  is 

accepted.  If  the  energy  increases  by  less  than  an  amount 
proportional  to  exp(-E/kT),  where  E  is  the  energy,  k  is 
Boltzmann's  constant  and  T  is  temperature,  the  move  is 
accepted.  Otherwise  the  move  is  rejected  and  the  system 
returned  to  its  previous  configuration.  If  a  statistically 
significant  number  of  moves  are  used,  average  values  of 
state  variables  calculated  by  MC  are  equivalent  to  those 
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found  by  molecular  dynamics.  However,  since  the  equations 
of  motion  are  not  solved,  monte  carlo  techniques  are  used 
only  to  find  static  properties.  (Haile,  1980:6-7) 

DYNAMO .  DYNAMO  is  the  atomistic  simulation  computer 
code  used  for  this  study.  DYNAMO  uses  molecular  statics  to 
perform  energy  minimization  calculations  (molecular 
statics),  and  molecular  dynamics  or  monte  carlo  (MC)  methods 
to  study  atomic  motion.  The  user  specifies  parameters,  such 
as  desired  temperature,  length  of  run  and  operating  method 
(.MS,  MC  or  MD),  in  an  input  file  read  by  DYNAMO.  DYNAMO 
uses  a  fifth-order  pred i ctor-corrector  routine  to  solve  the 
equations  of  motion  in  molecular  dynamics  numerically. 

Since  the  MD  routine  steps  forward  a  finite  time  in  each 
step,  the  behavior  of  the  system  and  its  properties  (energy, 
volume,  etc.)  with  respect  to  time  can  be  studied.  During 
each  step  of  a  monte  carlo  run,  a  random  atom  is  selected 
and  moved  a  certain  distance  which  depends  on  the  input 
temperature  of  the  system.  After  the  atom  is  moved,  the 
total  energy  of  the  system  is  calculated.  If  the  energy  of 
the  system  decreases,  or  increases  within  specified 
tolerances,  the  move  is  allowed;  otherwise  the  position  of 
the  atom  remains  unchanged.  Monte  carlo  is  used  to 
calculate  the  properties  of  canonical  ensembles. 

The  user  may  also  determine  if  a  new  lattice  is  created 
for  the  run  by  specifying  the  unit  cell,  atom  positions 
within  the  cell,  type  of  atom  at  each  position  in  the  cell 
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and  number  of  unit  cells  in  each  direction.  If  continued 
study  of  previous  work  is  desired,  a  restart  file  from  a 
previous  run  may  be  read  by  the  code  and  used  as  the  initial 
atomic  configuration,  instead  of  creating  a  lattice  for  a 
new  run.  Restart  files  are  include  in  the  standard  output 
from  a  DYNAMO  run  and  contain  information  about  the  system 
at  the  end  of  the  run  in  which  the  file  is  created.  The 
total  r. umbei  of  atoms,  atom  types,  physical  size  of  the 
system  and  the  positions  and  velocity  of  each  atom  are 
contained  in  the  restart  file. 

Another  output  from  DYNAMO  which  is  extremely  important 
to  this  work,  is  the  trajectory  file.  It  contains  the  same 
data  stored  in  the  restart  file,  but  the  trajectory  file 
saves  this  information  at  several  times  during  the 
simulation.  By  plotting  the  trajectories  of  an  atom  over 
several  thousand  time  steps,  one  may  make  qualitative 
judgements  about  the  phase  of  the  material.  That  is,  a 
liquid  will  move  randomly  and  over  relatively  great 
distances  whereas  a  solid  will  remain  within  a  short 
distance  of  its  original  position. 

The  third  output  from  DYNAMO  is  the  printout  file. 

This  file  stores  the  properties  of  the  system  at  specified 
intervals.  The  properties  include  the  time,  total  system 
energy  during  the  time  step,  and  the  instantaneous  and 
average  values  of  kinetic  energy,  potential  energy, 
temperature,  volume  and  pressure. 
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The  crystal  used  in  this  simulation  is  an  intermetal  1 ic 


compound,  CuT i .  The  CuT i  is  a  body-centered-tetragonal 
crystal,  which  can  be  visualized  as  a  cube  with  one  atom 
centered  inside  a  box  and  eight  other  atoms  placed  so  their 
centers  lie  on  the  corners  of  the  box.  For  this  simulation, 
individual  layers  of  atoms  lie  in  the  x-y  plane  and  are 
stacked  in  the  z-direction,  two  copper  planes  lying  on  2 
titanium  planes  (see  Figure  1).  The  choice  of  crystal 
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Figure  1,  CuTi  Crystal 


size  is  not  arbitrary;  the  system  needs  to  be  large  enough 
to  study  the  interface,  yet  small  enough  to  minimize 
required  computer  time.  A  long,  thin  crystal  is  chosen  to 
reduce  boundary  effects  along  the  x-axis  while  observing  the 
interface.  Y-  and  z-axis  boundary  effects  are  reduced,  but 
not  necessarily  eliminated,  by  making  the  crystal  periodic 
in  the  y  and  z  directions.  A  system  consisting  of  1080 
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atoms  (30  atomic  layers,  each  6  by  6  atoms)  is  studied  since 
the  code  would  not  need  to  be  modified  to  accommodate  this 
size  crystal  and  it  is  long  enough  to  allow  study  of  the 

i 

interface  over  meaningful  times  (>20000  steps).  If  the 
interface  crosses  the  crystal  too  quickly  little  or  no 
useful  data  is  obtainable. 

i 

The  Cu-Ti  system  is  created  using  the  lattice  creation 
capability  of  DYNAMO.  The  unit  cell  consists  of  four 
titanium  atoms  at  the  bottom  "corners"  of  the  cube,  a  copper 

► 

atom  centered  among  the  four  Ti  atoms  in  the  x-y  plane  and 
one  atomic  layer  above  them  in  the  z-direction  the  titanium. 
Two  similar  layers  with  the  atom  types  reversed  are  placed 
above  the  copper  atom  to  complete  the  unit  cell.  The  system 
used  is  formed  by  forming  a  block  of  30x6x3  unit  cells  of 
CuTi  . 

i 

Analysis  Methods 

Both  qualitative  and  quantitative  techniques  were  used 
to  analyze  the  output  from  the  simulation  runs.  Qualitative 
methods,  including  trajectory  plots  and  atom  position  plots, 
were  used  to  make  a  rough  estimation  of  the  interface 
location  and  the  chemical  order  or  disorder  of  the 
recrystallized  CuTi.  Mean  square  displacement,  atomic 
density  distributions  and  order  parameters  were  used  to 
quantify  the  position  of  the  interface  and  the  phase  of  the 
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material  in  the  melt  area.  A  comparison  of  the  analysis 
methods  (excluding  the  mean  square  displacement)  is  in 
Appendix  A. 

Trajectory  Plots.  Trajectory  plots  use  the  data  in 
trajectory  files  to  produce  a  picture  of  the  path  an  atom 
follows  during  the  run.  The  user  can  specify  the  time 
during  which  the  atoms  are  followed,  the  viewing  direction 
and  a  slice  of  the  system  to  view. 


Figure  2,  Sample  Trajectory  Plot 

A  sample  trajectory  plot  is  shown  in  Figure  2.  The  dots  on 
the  left  hand  side  of  Figure  2  are  atoms  held  in  place  by 
the  code.  The  larger  spots  to  the  right  of  the  fixed  atoms 
are  atoms  in  a  solid  phase.  Their  small  trajectories  result 
from  vibrational  motion.  At  higher  temperatures  the  spots 
will  get  larger,  at  lower  temperatures  the  spots  will 
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decrease  in  size.  The  area  right  of  the  solid  is  liquid, 
evidenced  by  the  long  pathlengths  and  lack  of  visible 
structure . 

Trajectory  plots  are  the  easiest  of  the  analysis 
methods  to  read,  but  they  must  be  used  with  caution.  If  the 
time  span  used  in  the  plot  is  too  small  or  the  time  steps 
too  large,  atoms  that  are  in  liquid  phase  may  appear  to  be 
an  amorphous  solid  because  there  are  not  enough  points  to 
plot.  The  atom  in  the  liquid  will  seem  to  move  only  a  short 
distance.  On  the  other  hand,  viewing  the  trajectories  over 
too  much  time  is  equally  undesirable.  Similar  care  must  be 
taken  when  selecting  the  slice  thickness.  Too  thin  a  slice 
will  result  in  a  sparse  plot  revealing  little  information. 

If  the  slice  is  too  thick,  the  overlapping  trajectories  of 
atoms  along  the  viewing  axis  will  smear  the  definition  of 
the  plot.  Because  of  these  limitations,  trajectory  plots 
are  only  used  for  a  crude  estimate  of  interface  position. 

Atom  Position  Plots.  Atom  position  plots,  also  called 
atom  projection  plots,  show  a  snapshot  image  of  the  system. 
These  plots  use  the  atom  position  information  contained  in 
the  restart  files  to  plot  the  atoms  in  their  relative 
positions.  The  user  may  specify  the  viewing  angle,  the 
orientation  of  the  system,  and  the  depth  and  thickness  of 
the  slice  viewed.  A  sample  projection  file  is  shown  in 
Figure  3. 
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Figure  3,  Sample  Projection  Plot 

The  filled  circles  in  the  projection  plot  represent  titanium 
atoms,  the  open  circles  are  copper.  Moving  from  the  left 
boundary  to  the  right,  three  distinct  regions  are  seen  in 
the  crystal.  The  three  layers  on  the  left  border  are  held 
in  place  by  the  code  as  shown  by  the  rigid  structure.  Order 
is  evident  in  the  next  six  atomic  planes  to  the  right,  but 
much  looser  because  the  atoms  have  moved.  Antisite  defects, 
when  two  atoms  of  different  kind  change  places,  are  also 
present,  an  example  being  the  fifth  atomic  layer  from  the 
left.  At  the  top  of  this  layer,  a  copper-titanium  pair  are 
seen  where  a  Ti-Ti  pair  are  expected. 

Projection  plots  are  useful  to  see  defects  in  the 
crystal,  such  as  vacancies  or  antisite  defects,  and  to  see 
if  the  system  shows  structure.  The  plots  are  not  useful  to 
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determine  the  position  or  movement  of  a  liquid-solid 
interface  because  it  is  not  possible  to  determine  if  a 
disordered  region  is  a  liquid  or  an  amorphous  solid. 

Mean  Square  Displacement  (MSP).  An  atom’s  mean  square 
displacement  is  the  square  of  the  vector  distance  the  atom 
moves  from  its  starting  position  over  a  specified  time.  To 
calculate  the  mean  square  displacement,  the  system  is 
divided  into  groups  of  atoms,  each  one  atomic  layer  or  1.5 

o 

angstroms  ( A )  thick,  spanning  the  y-z  plane.  The  user 
specifies  the  trajectory  file  to  be  used  and  the  starting 
and  ending  times  of  interest.  The  squared  displacements  of 
all  atoms  in  the  group  are  summed  and  divided  by  the  number 
of  atoms  in  the  group  at  the  initial  time.  If  an  atom  from 
group  12  crosses  the  boundary  into  group  13  during  the  MSD 
calculation,  it  is  still  considered  to  be  part  of  group  12 
during  the  remainder  of  the  calculation. 

The  MSD  is  a  method  of  determining  the  position  of  the 
interface  and  the  phase  of  the  system.  The  self-diffusion 
coefficient  for  liquids  is  several  orders  of  magnitude 
higher  than  for  solids,  so  atoms  in  liquids  move  farther 
than  in  solids.  In  this  study,  the  mean  square  displacement 
is  averaged  over  several  groups  in  the  apparent  liquid 
region  in  a  trajectory  plot.  The  location  of  the  interface 
is  defined  as  the  position  where  the  MSD  is  one-half  the 
average  group  mean  square  displacement  of  the  liquid.  To 
identify  the  phase  of  the  region,  averagf  group  MSDs  are 
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calculated  for  the  groups  in  the  region  as  a  function  of 
time.  The  MSD  will  approach  a  constant  value  over  time  for 
solids,  but  will  increase  linearly  with  time  for  liquids 
(Hai 1 e,  1988 ,  36) . 

Atomic  Density  Distribution.  Atomic  density 
distributions  show  the  distribution  of  atoms  along  the 
x-axis  of  the  system.  The  system  is  broken  into  smaller 
groups  than  for  a  mean  square  displacement,  each  having  a 

o 

thickness  of  .375  A.  The  number  of  atoms  in  each  group  is 
plotted  against  the  midpoint  of  the  group.  If  the  system  is 
crystalline,  the  distribution  will  have  high,  sharp  peaks 
one  or  two  groups  wide,  and  deep  valleys  of  very  few  to  no 
atoms.  An  amorphous  or  liquid  system  will  be  more  level 
than  the  crystal . 


Order  Parameter.  The  order  parameter  is  a  way  of 
measuring  the  crystalline  nature  of  a  system.  The  order 
parameter  is  found  by 


P  gk 


1 

AtN(N  -  1 ) 


(12) 


p9„  is  the  order  parameter  of  group  g  in  the  k  direction,  N 

is  the  number  of  atoms  in  group  g,  \k  is  the  distance 
between  atomic  layers  in  the  k  direction  (the  lattice 
constant),  and  r£-r(  is  the  distance  from  atom  i  to  j  in  the 
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k  direction.  At  is  the  number  of  trajectory  files  used  when 
calculating  an  average  order  parameter  for  a  given  time 

during  a  run. 

i 

As  a  system  becomes  crystalline,  the  average  distance 
between  atoms  in  a  given  direction  will  approach  the  lattice 
constant  for  that  direction,  \k.  For  a  perfect  crystal  p*-l 

> 

for  all  directions.  if  the  crystal  is  heated,  the  order 
parameter  will  be  slightly  less  than  one,  since  the  atoms 
will  move.  The  interface  between  crystalline  and  amorphous 
regions  is  defined  here  as  <pyi>-0.2,  <pyi>  is  the  average  of 
the  order  parameters  in  the  y-  and  z-directions. 
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Procedure 


Overview 

The  type  of  approach  used  was  required  not  only  because 
of  the  limited  research  time,  but  also  due  to  the  limited 
amount  of  time  available  on  the  Ohio  Supercomputing  Center 
Cray  Y-MP.  The  crystal  was  created  and  heated  to  the 
temperatures  to  be  used  in  the  freezing  runs.  The 
liquid-solid  interface  was  created  using  a  modified  monte 
carlo  subroutine  and  the  system  was  cooled  various 
temperatures.  The  output  was  analyzed  and  decisions  to 
either  continue  at  the  current  temperature  or  to  lower  the 
temperature  were  made. 

Liquid-Solid  I nterf ace  Creation 

After  the  crystal  was  created,  it  was  heated  to  1000  K, 
1100  K  and  1200  K  to  prepare  it  for  the  next  phase:  melting 
part  of  the  crystal  to  set  up  the  interface.  The  system  was 
minimized  after  heating  to  restore  the  Cu-Ti  to  a  perfect 
crystal  lattice.  The  y-  and  z-boundaries  were  fixed  at 
their  post-minimization  values  so  no  shear  stresses  will 
occur  between  the  held  and  free  atoms  when  the  atoms  are 
released  for  freezing  runs.  The  x-boundary  was  extended  to 

o  O 

90  A,  creating  a  40  A  void  along  the  x-axis  to  allow  the 
system  to  melt.  If  the  x-boundary  was  fixed  at  its 


minimized  value,  the  system  would  not  melt  because  the 
system  could  not  expand,  thus  maintaining  the  atoms  in  their 
lattice.  The  crystal  was  now  ready  to  create  the  interface. 

The  liquid-solid  interface  was  made  using  the  monte 
carlo  routine  in  DYNAMO,  modified  as  follows.  Atoms  with  an 

o 

x-coordinate  less  than  9.0  yi  were  held  in  place.  Melting 
was  accomplished  at  3000  K,  much  higher  than  the 
exper imental  melting  temperature,  to  reduce  the  time 
required  to  set  up  the  interface.  The  output  from  the 
melting  run  is  hereafter  referred  to  as  "the  melt". 

Resolidification  Using  Mo  1  ecu  1  ar  Dynami c.s 

After  the  liquid-solid  interface  was  set  up,  molecular 
dynamics  was  used  to  resolidify  the  system.  The  boundaries 
remain  fixed  at  the  same  values  as  the  melting  runs.  The 
first  three  atomic  layers  from  the  x-axis  are  held  to 
simulate  a  solid  CuT i  crystal  extending  in  the  negative 
x-direction.  If  the  atoms  are  not  held,  the  system  may  melt 
completely  because  it  sees  a  void  for  x  less  than  zero  due 
periodicity  and  the  large  void  introduced  in  the  melting 
procedure  (total  system  melting  was  observed  in  the  initial 
recrystal  1 izat ion  attempts). 

The  first  step  in  recrystallization  was  to  reduce  the 
temperature  to  the  value  corresponding  to  the  boundary 
lengths  used  during  melting.  This  was  accomplished  using 
only  2000-4000  steps  to  lessen  recrystallization  during  the 


cooling.  In  all  the  runs  described  here,  a  time  step  was 
defined  as  2xl0~ls  seconds.  After  the  crystal  was  cooled, 
several  runs  are  completed  over  longer  times,  10000  or  20000 
steps.  After  each  run  was  completed,  the  output  was 
examined  to  determine  the  location  of  the  interface,  the 
order  of  the  atoms  in  the  solid,  and  the  phase  of  the 
disordered  volume  (liquid  or  amorphous  solid).  If  the 
interface  has  not  moved  noticeably,  the  temperature  was 
lowered  100  K,  and  the  system  run  again.  The  means  of 
analysis  are  covered  in  the  following  chapter. 

Recrvstal 1 izat ion  Attempts.  Based  on  the  work  by  Wood 
and  the  experimental  value  of  the  melting  point  of  CuT i , 

1200  K  was  the  temperature  used  for  the  first  solidification 
attempt.  The  crystal  was  held  at  1200  K  for  20000  time 
steps.  The  subsequent  plot  of  the  trajectories  showed  the 
crystal  had  completely  melted.  A  similar  set  of  runs  made 
at  1100  K  also  produced  a  melted  system,  so  a  decision  was 
made  to  hold  three  atomic  layers  during  the  crystallization 
runs.  The  rationale  for  immobilizing  the  layers  was  that 
the  interface  surface  under  study  may  be  of  such  a  high 
energy  state  that  the  system  found  it  easier  to  melt  than  to 
maintain  and  build  on  the  interface.  Since  monte  carlo  was 
not  used  for  recrystallization,  the  forces  on  each  atom  and 
the  atoms'  velocities  were  set  to  0  during  the  force 
calculations  in  molecular  dynamics. 
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Starting  from  the  melt  at  3000  K,  the  1200  K 
recrystallization  attempt  was  repeated.  The  system  was 
cooled  to  1200  K  in  two  separate  runs:  the  first  of  1000 
steps  to  2000  K  and  the  second  to  1200  K  in  1500  steps. 

After  cooling,  the  CuT i  was  held  at  1200  K  for  2  subsequent 
runs,  one  10000  steps  and  one  5000  steps.  Examination  of 
the  atom  trajectories  showed  no  noticeable  movement  of  the 
liquid-solid  interface.  A  plot  of  the  atoms  showed 
structure  in  the  three  rows  held  and  in  the  two  adjacent 
rows.  The  structure  in  the  two  free  rows  is  seen  because  of 
the  binding  effect  of  the  stationary  rows.  The  structured 
region  is  limited  to  two  rows  due  to  the  effective  range  of 
the  interatomic  potential.  if  the  system  would 
recrystallize  at  1200  K  it  would  take  too  long  to  see  in 
terms  of  simulation,  and  hence  computer,  time.  The  time  of 
the  simulation  at  1200  K  was  35  picoseconds  (ps). 

The  temperature  was  lowered  to  1100  K  in  an  attempt  to 
recrystallize  the  system  faster.  In  addition  to  the  lower 
temperatures,  the  boundary  lengths  in  the  y-  and 
z-directions  were  shortened  to  the  edge  lengths  of  an 
equilibrated  1100  K  crystal.  Again  the  runs  were  started 
from  the  melt  at  3000  K  and  reduced  to  1100  K  in  2  runs  of 
2000  steps  each.  The  initial  cooling  runs  were  longer  to 
allow  the  shock  introduced  by  squeezing  the  system  to  pass 
through  the  system.  Three  10000  step  computer  runs  were 
performed  on  the  system  at  1100  K.  Analysis  of  the 


trajectories  and  atom  positions  failed  to  show  quantifiable 
movement  of  the  interface,  even  though  the  system  had  been 
run  for  68  ps. 

The  temperature  was  again  lowered,  to  1000  K,  and  the 
boundaries  reduced  accordingly.  To  reduce  the  computer  time 
required  in  the  initial  cooling  phase,  the  restart  file  from 
the  first  10000  step  run  at  1100  K  was  chosen  as  a  starting 
configuration  for  the  1000  K  runs.  No  initial  cool-down 
runs  were  performed  on  the  CuT i ;  instead,  the  first  run  at 
the  new  temperature  lasted  10000  steps.  A  quick  overview  of 
the  trajectory  plots  revealed  a  small  area  near  the  z  =  0 
plane  appeared  to  crystallize  early  in  the  run,  between  2000 
and  6000  timesteps,  only  to  disappear  by  the  end  of  the  run. 

he  crystal  appeared  to  form  facets  into  the  liquid  instead 
of  growing  perpendicular  to  the  interface  surface  (see 
Figure  4).  This  is  similar  to  the  faceting  observed  in 
silicon  by  Landman  in  the  simulation  of  Si  (111)  (Landman, 
et  al . ,  1988b).  Although  the  crystal  seemed  to  melt  by  the 
end  of  the  run,  its  appearance  was  encouraging  enough  to 
continue  at  1000K  for  20000  more  steps.  The  additional  time 
at  T  =  1000  K  did  not  produce  the  desired  results.  More 
faceting  was  observed,  but  the  facets  would  melt  before  they 
had  grown  large  enough  to  sustain  recrystallization.  More 
runs  at  1000  K  were  deemed  unlikely  to  recrystallize  the 


system. 


4^: 

.:*v 


mm 


Figure  4,  Trajectories,  T  =  1000  K 
t  =  32  -  36  ps 


The  next  temperature  tried  was  900  K,  with 
corresponding  boundary  lengths.  The  CuTi  was  held  at  900  K 
for  20000  timesteps  using  the  output  from  the  first  1000  K 
run  as  a  starting  file.  At  the  end  of  this  run  (88  ps  from 
the  initial  melt),  two  atomic  layers  of  crystal  appeared  to 
have  grown  and  remained  at  the  end  of  the  run.  The  system 
was  held  at  900  K  for  129  ps  more,  which  formed  only  two 
additional  atomic  layers  of  crystal.  More  crystal  growth 
had  been  expected.  The  temperature  was  lowered  100  K  once 
more  in  an  attempt  to  freeze  the  system  in  an  amorphous 
state. 

The  CuTi  was  held  at  800  K  for  40  ps.  The  starting 
point  was  the  end  of  the  first  1000  K  run,  48  ps  from  the 
melt.  After  a  the  additional  40  ps,  the  projection  plot  of 
the  final  configuration  showed  no  obvious  movement  of  the 
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interface  and  trajectory  plots  indicated  that  the  atoms  in 
the  melted  region  seemed  to  move  less  as  the  simulation 
progressed.  A  plot  of  MSD  vs  time  for  the  disordered  region 
(Figure  5)  showed  a  linear  increase  in  the  MSD.  It  is  not 
evident  whether  this  is  an  amorphous  solid,  glass  or  liquid, 
since  glasses  do  exhibit  some  liquid-like  properties. 

After  failing  to  produce  an  amorphous  solid  at  800  K, 
the  three  layers  of  atoms  that  had  been  held  were  released. 
The  output  from  the  final  900  K  run,  217  ps  from  the  melt, 
was  held  for  an  additional  40  ps.  Although  S00  K  was  the 
intended  temperature,  the  system  was  actually  run  at  818  K. 
Projection  plots  and  trajectory  plots  still  did  not  show  the 
interface  moving.  Comparing  the  order  parameters  at  8  ps 
intervals  through  the  run  confirmed  the  interface  had  not 
noticeably  moved. 


Figure  5,  MSD  vs  Time,  T  =  800  K 


A  run  at  773  K  was  set  up  and  run  for  80  ps  timesteps 
in  a  final  attempt  to  form  an  amorphous  solid.  However, 
analysis  of  the  trajectory  plots  and  projection  plots  showed 
a  definite  crystal  growth  in  the  system.  The  CuTi  was  run 
for  an  additional  80  ps  at  773  K,  and  further  crystal  growth 
was  noted.  Figures  6a-c  show  the  projection  plots  of  the 
system  at  the  beginning  of  the  runs  at  773  K,  after  80  ps  at 
773  K  and  after  160  ps,  respectively. 


Figure  6a,  Atom  Positions  at 
Beginning  of  773  K  Runs 


O  O  O  • 

@  ®  <Q  ' 

0  $  C  0® 

<d  o  e 

O  <D 

©  © 

•rP.O  <S> 


-  n  oq._  9  X® 


.o. 


O  o 

<0  o  O  Q  *  & 

•  •  ® 

O  Q)  0  O  C 


5<*%> 


^T\ 


•  c?  •  cb>oo  »00^“d 


o 
o 


Figure  6b,  Atom  Positions  after 
80  Picoseconds  at  773  K 
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Figure  6c,  Atom  Positions  after 
160  Picoseconds  at  773  K 


The  crystal  clearly  grew  from  left  to  right,  starting  about 
one-quarter  of  the  distance  across  the  crystal  and  ending 
almost  two-thirds  across.  The  trajectory  plots 
corresponding  to  each  projection  plot  are  shown  in  Figures 
7a-c . 
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Figure  7a,  Trajectory  Plot 
0-16  Picoseconds,  T  =  773  K 
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Figure  7b,  Trajectory  Plot 
80-86  Picoseconds,  T  =  773  K 
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Figure  7c,  Trajectory  Plot 
144-160  Picoseconds,  T  =  773  K 


The  trajectory  plots  also  show  the  crystal  grew  from  left  to 
right,  covering  almost  half  the  crystal  in  160  picoseconds. 
The  output  from  these  runs  was  analyzed  in  greater  detail  to 
determine  the  interface  velocity.  This  analysis  is  covered 
in  Chapter  IV. 
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As  seen  in  Figures  6  and  7,  the  CuT i  appears  to  have 

recrystallized  across  a  large  portion  of  the  system.  To  get 

more  precise  values  for  the  extent  of  recrystallization  and 

the  velocity  of  the  interface,  more  detailed  analysis  is 

required  using  atomic  density  plots  and  the  order  parameter. 

MSD  calculations  were  unreliable  for  the  818  K  and  773  K 

runs.  Before  this  analysis  was  done,  an  approximation  of 

the  interface  velocity  was  made  using  the  MSD  output  from 

the  initial  1000  K  run  and  the  818  K  run.  An  average  mean 

square  displacement  was  calculated  over  ten  groups  in  the 

liquid  region  for  five  time  spans.  The  diffusion 

coefficient  was  calculated  for  each  At  using  Equation  11, 

and  an  average  diffusion  coefficient  was  found  for  the 

liquid  region.  The  interface  velocity  was  calculated  using 

0 

Equation  10  using  X-1.5/!.  The  results  are  shown  in  Tables 
1  and  2.  The  velocity  limit  at  818  K  is  lower  than  the 
limit  at  1000  K.  This  may  be  due  to  reasons  covered  in  the 
chapter  on  previous  work. 


Table  1 

Interface  Velocity  Limit,  1000  K 


|  At  (ps) 

a  2 

MSD  (A  ) 

D  (cm  2 / s ) 

Interface  1 

Velocity  I 

Limit  (m/s)  f 

1  8 

|  12 

1  16 

I  20 

|  Average 

4.7±1 .2 
8.9±1  .2 

1 2 . 0+2 . 2 

1 3 . 3±2 . 2 

20 . 7+3 . 0 

2 . 0± .  5  x  10-5 
1.9+.4  x  10-5 

1 . 7±.  3  x  10-5 
1.4+.3  x  10-5 

1 .7±.3  x  10-5 

1 . 7±. 3  x  10-5 

11.3+2.0  | 

Table  2 

Interface  Velocity  Limit,  818  K 


B  At  (ps) 

0  2 

MSD  ( A  ) 

D  (cm  2/s) 

1 nterf ace 
Velocity 
Limit  (m/s) 

8 

16 

24 

32 

40 

Average 

3 . 9± 1 . 1 
7.3+1 .4 

1 1 .5±1 .9 

1 7 . 7  +  2 . 9 

1 8 . 8±4 . 1 

8. 2  +  2. 3  x  10-6 
7.6+1 .6  x  10-6 

8 . 0± 1 . 3  x  10-6 
9.2+1 . 5  x  10-6 

7 . 8± 1 .7  x  10-6 
8. 2+1 .7  x  10-6 

5.4±1  .  1 

Trajectory  and  Projection  Plots.  773  K 

Examining  Figures  6  and  7,  it  appears  the  interface  has 
moved  approximately  five  atomic  layers  or  7.5  angstroms  to 
the  right  during  the  runs.  A  closer  examination  of  the 
projection  plots  shows  that  as  the  system  recrystallizes  it 
does  not  grow  back  to  its  initial  configuration.  Several 
anti-site  defects  are  visible,  but  no  point  defects  appear. 


The  anti-site  defects  occur  because  atoms  entering  incorrect 
lattice  sites  do  not  possess  enough  energy  to  break  free, 
thus  allowing  a  correct  type  of  atom  to  enter.  The 
trajectory  plots  in  Figure  7  show  small  trajectories  for  the 
atoms  in  the  disordered,  as  well  as  the  ordered  region. 

This  could  have  two  causes:  1)  the  time  span  of  the 
trajectory  plots  was  too  short,  or  2)  the  disordered  region 
has  glass-like  properties.  If  2)  is  correct,  the  high 
temperature  would  cause  the  system  to  be  viscous  like  molten 
glass  and  still  allow  diffusive  crystal  growth.  The 
trajectory  and  projection  plots  can  only  give  a  rough  idea 
of  the  growth  of  the  crystal;  the  other  analysis  techniques 
were  used  to  calculate  the  interface  position  and  velocity. 

Atomic  Density  Plc.ts 

The  number  of  atoms  within  a  given  group  was  plotted 
versus  the  center  x-coordinate  of  the  group  at  the  beginning 
and  end  of  the  runs  at  773  K.  The  plots  are  shown  in 


Figures  8a-b. 
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Figure  8a,  Atomic  Density 
t  =  0  ps,  T  =  773  K 


Figure  8b,  Atomic  Density 
t  =  160  ps,  T  =  773  K 


As  the  interface  moves  to  the  right,  the  number  of  sharp 
peaks  and  deep  valleys  increases.  Each  peak  roughly 
corresponds  to  the  position  of  a  crystal  layer.  As  the  plot 
moves  to  more  disordered  regions,  the  difference  between  the 
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number  of  atoms  in  the  peaks  and  valleys  becomes  less.  The 
atomic  density  plots  confirm  that  the  system  is  growing  in 
layers  and  not  forming  an  amorphous  solid.  The  location  of 
the  interface  is  found  using  the  order  parameter. 

Interface  Location  Using  Order  Parameter 

The  order  parameter,  p,  was  averaged  for  every  eight 
picosecond  span  for  the  two  773  K  runs.  That  is,  the 
average  order  parameters  between  0  and  8  ps,  8  and  16  ps,  16 
and  24  ps  and  so  on,  were  calculated.  The  order  parameter 
was  calculated  using  Equation  13  for  thirty  five  groups  with 
widths  of  1.5  angstroms.  For  analysis  purposes,  it  is 
assumed  the  average  order  parameter  occurred  at  the  midpoint 
of  each  time  span  and  each  group.  The  y  and  z  order 
parameters  were  averaged  for  groups  bordering  the  interface. 
The  location  of  the  interface  was  calculated  using  linear 
interpolation  to  find  <pyi>  =  0.2.  The  interface  was 
tracked  for  the  duration  of  both  runs  and  plotted  vs  time  in 
Figure  9.  Figure  9  shows  the  interface  does  not  move  with  a 
constant  velocity,  similar  to  what  Landman  et  al .  observed 
for  the  recrystallization  along  the  flat  surface  of  silicon. 
The  interface  stays  in  the  same  vicinity  for  the  first  52  ps 
grows  at  a  very  fast  rate  in  the  next  32  ps,  remains 
constant  for  a  short  period  and  the  grows  again.  The 
periods  where  no  growth  occurs  may  occur  after  one  or  two 


layers  grow  then  the  atoms  within  those  layers  rearrange 
themselves  to  a  lower  energy  configuration.  The  position  of 
the  interface  remains  constant  until  the  crystal  has  stopped 
this  rearrangement . 
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Figure  9,  Interface  Position  vs  Time 
T  =  773  K 


The  positions  at  specific  times  during  the  run  are  used 
to  calculate  the  velocity  of  the  interface  between  two 
times.  Time  t  =  0  is  the  beginning  of  the  first  run  at  773 
K.  The  velocity  is  calculated  simply  by  dividing  the 
distance  the  interface  moved  (9.98  angstroms)  by  the  time 
(152  ps).  The  average  velocity  of  the  interface  over  152 
picoseconds  is  6.6  m/s,  on  the  order  expected  in 
diffusion-limited  growth.  During  the  periods  of  rapid 
movement  the  interface  reached  speeds  of  up  to  18  m/s 
(between  124  and  152  ps),  but  longer  runs  are  needed  to  see 


if  these  speeds  are  sustained  for  long  times. 

The  simulated  interface  speed  is  .02  m/s  above  the 
velocity  limit  calculated  earlier  for  818  K.  This  actually 
agrees  very  well  with  theory,  because  the  values  for  K  are 
not  well  defined.  The  system  could  be  recrystallizing 
faster  than  expected  since  the  atoms  do  not  have  enough 
energy  to  diffuse  out  of  incorrect  lattice  sites  and  do  not 
move  before  another  layer  freezes  over  them.  The  interface 
would  not  halt  for  long  periods  while  the  atoms  in  the 
crystal  reorganized  to  a  lower  energy  configuration. 

What  is  interesting  to  note  is  the  system  refreezes  in 
a  chemically  disordered  state.  If  the  low  speeds  are  the 
result  of  diffusion-limited  growth,  the  crystal  should 
refreeze  into  its  original  configuration,  since  the 
recrystallization  is  slowed  as  atoms  diffuse  into  the  proper 
locations.  If  the  type  of  atom  entering  the  lattice  site 
does  not  make  any  difference,  as  is  apparently  seen  in  the 
simulation,  then  the  interface  velocity  should  be  on  the 
order  of  col  1 ision- 1 imited  growth  (~  100  m/s).  If  an  atomic 
layer  requires  a  certain  proportion  of  atom  types  be  present 
in  that  layer  but  not  In  a  specific  configuration,  a 
diffusive  type  growth  would  be  observed.  No  analysis  was 
done  to  determine  if  this  is  what  occurred. 
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Conclusions 


The  recrystallization  and  behavior  of  the  liquid-solid 
interface  of  a  CuT i  body-centered-tetragonal  crystal  were 
studied  using  atomistic  simulation.  An  interface  was  set  up 
in  the  crystal  using  monte  carlo  and  the  system  was  cooled 
at  various  temperatures  using  molecular  dynamics. 

Chemically  disordered  recrystal  1 i zat ion  occurred  at  773  K: 
sustained  recrystallization  was  not  observed  at  higher 
temperatures  and  lower  temperatures  were  not  simulated. 

Recrystallization  appeared  to  occur  as 
diffusion- 1 imited  growth  at  773K  based  on  the  magnitude  of 
the  interface  velocity  taken  from  the  simulation.  If  the 
growth  was  col  1 ision- 1 imited,  the  speeds  would  be  an  order 
of  magnitude  higher.  The  recrystallization  at  this 
temperature  is  not  truly  diffusion-limited  growth  or  the 
number  of  antisite  defects  would  be  significantly  reduced, 
possibly  to  zero. 

Although  the  growth  appeared  to  be  diffusion-limited, 
the  regrown  crystal  was  observed  to  be  chemically 
disordered.  The  antisite  defects  arise  due  to  the  lack  of 
mobility  of  the  atoms  in  the  amorphous  region  at  the  low 
temperature.  Attempts  to  recrystallize  the  CuT i  in  its 
original  configuration  must  be  made  at  higher  temperatures 


40 


to  give  the  atoms  in  the  liquid  phase  the  energy  required  to 
move  out  of  -improper  lattice  sites  Lefcre  other  layers  form 
over  them. 

By  no  means  was  this  an  exhaustive  study  of  the 
behavior  and  properties  of  the  liquid-solid  interface  and 
recrystallization  of  CuT  i ,  primarily  because  of  time 
constraints.  The  lack  of  previous  work  on  metal  alloys  was 
also  hindering  since  no  obvious  starting  point  was 
available.  Using  molecular  dynamics  to  simulate 
recrystallization  of  a  concentrated  metal  alloy  like  CuT i 
was  meaningful  since  no  work  has  been  published  on  it.  The 
velocities  resulting  from  the  simulation  at  773  K  are  within 
20%  of  the  limit  at  818  K,  which  is  good  agreement 
considering  the  lattice  constant  is  not  well  defined. 

Several  questions  arose  during  the  study,  such  as:  1)  If 
the  recrystal  1 ization  is  governed  by  dif fusion- 1 imi ted 
growth,  why  are  so  many  anti-site  defects  seen?  2)  What  is 
happening  during  the  periods  when  the  interface  does  not 
move?  or  3)  What  is  the  thickness  of  the  interface?  A  lack 
of  both  computer  time  and  real  time  did  not  allow  deeper 
investigation  of  these  questions,  but  this  study  does 
provide  a  foundation  for  future  work  in  the  behavior  of 
liquid-solid  interfaces  of  CuTi. 


The  four  methods  used  to  determine  the  location  and 
movement  of  the  liquid-solid  interface  during  the  773  K  runs 
are  compared  here  to  show  they  agree  with  each  other. 

Figure  A-l  is  the  plot  of  atom  positions  after  80  ps  at  773 
K.  The  atoms  have  some  order  in  the  first  11  to  12  atomic 
planes  from  the  left  boundary,  beyond  the  12  plane,  no  order 

O 

is  evident.  If  one  atomic  plane  is  approximately  1.5  A 

o 

thick,  the  interface  is  between  16  and  18  A.  Figure  A-2  is 
a  trajectory  plot  of  the  system  between  80  and  96  ps  at  773 
K.  The  solid  region,  denoted  by  the  small  spots,  extends 

o 

about  11  atomic  layers,  or  about  17  A  into  the  system  from 
the  left.  A  plot  of  the  atom  densities  vs  distance  from  x  = 
0  is  seen  if  Figure  A-3.  The  sharp  peaks  and  valleys  on  the 
left  indicate  a  crystalline  region,  the  relatively  smooth 
region  is  either  amorphous  solid,  glass  or  liquid.  The 

O 

sharp,  high  peaks  end  approximately  at  x  =  17  A,  which 
suggests  that  is  the  position  of  the  interface.  The  final 
and  most  quantitative  way  of  showing  the  movement  of  the 
interface  is  the  order  parameter.  Table  A-1  shows  the  y  and 
z  order  parameters  after  80  picoseconds  at  773  K.  The  two 
order  parameters  are  averaged  and  the  interface  is  chosen 
where  the  average  is  0.2.  Although  the  average  order 
parameter  for  group  4  is  less  than  .2,  this  is  because  of 
movement  into  the  vacuum  to  the  left  of  the  system  not  the 
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interface.  The  interface  is  in  group  13,  or  between  17.5 
0 

and  19  A  However,  further  investigation  shows  the  center 

O 

of  mass  has  moved  approximately  1.8  A,  so  the  interface  is 

O 

actually  between  15.7  and  17.2  A. 
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Figure  A-1 

Atom  Positions  after  80  ps  at  773  K 
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Figure  A-2 

Atom  Trajectories,  t  =  80-96  ps 
T  =  773  K 
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Atomic  Density  at  t  =  80  ps 
T  =  773  K 


Table  A- 1 

Order  Parameters  after  80  ps,  T  =  773  K. 


Group 

Boundaries 

(A) 

y 

z 

1 

-0.5 

1.0 

0.26147 

0.11005 

2 

1.0 

2.5 

0.46975 

0.27542 

3 

2.5 

4.0 

0.48089 

0.34289 

4 

4.0 

5.5 

0.55569 

0.43370 

5 

5.5 

7.0 

0.45443 

0.43219 

6 

7.0 

8.5 

0.44540 

0.51101 

7 

8.5 

10.0 

0.51584 

0.45268 

8 

10.0 

11.5 

0.51741 

0.47067 

9 

11.5 

13.0 

0.46813 

0.44353 

10 

13.0 

14.5 

0.45740 

0.45532 

1 1 

14.5 

16.0 

0.38969 

0.40226 

12 

16.0 

17.5 

0.33424 

0.37536 

13 

17.5 

19.0 

0.20544 

0.18161 

14 

19.0 

20.5 

0.10207 

0.06441 

15 

20.5 

22.0 

0.00370 

0.01529 

16 

22.0 

23.5 

-0.00381 

-0.01328 

17 

23.5 

25.0 

-0.02029 

-0.00656 

T  ime 


Center  of  Mass 


x 


At  start  of  run 
After  80  ps 


25.73126 

27.54551 


y 

9.36897 

9.25278 


(A) 


z 

9.47101 

9.39467 
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